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ABSTRACT 



Min et al. (2009) presented two complementary techniques that use the diffusion approximation to allow efficient Monte-Carlo radi- 
ation transfer in very optically thick regions: a modified random walk and a partial diffusion approximation. In this note, I show that 
the calculations required for the modified random walk method can be significantly simplified. In particular, the diffusion coefficient 
and the mass absorption coefficients required for the modified random walk are in fact the same as the standard diffusion coefficient 
and the Planck mean mass absorption coefficient. 
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1. Introduction 

The problem of Monte-Carlo radiation transfer in very optically 
thick regions - such as in the midplane of circumstellar disks - 
is challenging. Without any approximations, photon packets can 
get trapped for millions of interactions, increasing th e required 
compu tational time by several orders of magnitude. [Min et ail 
(2009, hereafter M09) presented two complementary methods to 
greatly improve the efficiency of Monte-Carlo radiation transfer 
codes in very optically thick regions: a modified random walk 
(MRW) and a partial diffusion approximation (PDA). The MRW 
prevents photons from getting stuck in very optically thick re- 
gions, and the PDA allows temperatures to be calculated in re- 
gions that see few or no photons. 

The essence of the MRW method is that instead of comput- 
ing thousands to millions of individual absorption or scattering 
events for a single photon in these optically thick regions, one 
can make use of the solution to the diffusion approximation in- 
side small regions to propagate the photon efficiently. Monte- 
Carlo radiation transfer codes propagate photons in a grid made 
up of cells of constant density and temperature. Therefore if the 
mean optical depth to the edge of a cell is much larger than unity, 
one can set up a sphere whose radius is smaller than the distance 
to the closest wall, inside which the density will be constant, and 
travel to the edge of a sphere in a single step using the diffusion 
approximation. 

A probability distribution function is used to sample the true 
distance traveled to exit this sphere (since the photon would fol- 
low a random walk inside the sphere, rather than moving in a 
straight line). This true distance, which depends on the radius of 
the sphere and the local diffusion coefficient D, can then be used 
along with the mass absorption coefficient k to compute the to- 
tal amount of energy deposited in the dust during the diffusion. 
This is required in order to compute the temperature in the cell 
accurately. Finally, the photon is emitted from a random position 
on the surface of the sphere with a frequency sampled from the 
Planck function. 



M09 provide equations for D, k, and the dust emission coef- 
ficient 7] Y , taking into account that photons can be both scattered 
and absorbed and re-emitted. In M09, the suggested algorithm 
is to first calculate r\ v iteratively, and to then use rj v to compute 
D and k. In this note, I show that r\ v does not need to be solved 
iteratively, but can be solved directly, and I use this solution to 
show that D and k can in fact very easily be computed, result- 
ing in both a simpler implementation of the MRW, and in some 
cases performance gains. 

2. Derivation 

2. 1 . Emission coefficient 

In the presence of isotropic scattering, the emissivity of dust in 
local thermodynamic equilibrium (LTE) is given by 



T]y = KyBy(T) + CTyJy 



(1) 



where v is the frequency of the radiation, k v is the mass absorp- 
tion coefficient, B V (T) is the Planck function at the temperature 
T of the dust, o~ v is the mass scattering coefficient, and J Y is the 
mean intensity of the radiation field. Assuming that the radiation 
is isotropic, J v = 7 V , where 7 V is the intensity of the radiation. 
In the optically thick regime, I v - S v , where S v is the source 
function. Therefore, 



KyBy(T) + CTySy 



(2) 



The source function is defined as the ratio of the total emissivity 
to the total extinction, which in this case is 



TJy 

Xv 



(3) 



where Xv is the mass extinction coefficient (x v = « v + <x v ). 
Therefore, Equation (|2]i can be rewritten as 



TJy - Ky By(T) + CTy . 

Xv 



(4) 
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Re-arranging this equation, one obtains 

K V By{T) 

which can be simplified, since: 

Ky Ky Ky 



Therefore, 

T]y = XV By{T) 

The source function for this emissivity is 
S v = = Bv {T), 

Xv 



Xv 



(5) 



(6) 



(7) 



(8) 



which is expected for thermal emission from dust in LTE in the 
optically thick regime. 

2.2. Comparison with the M09 emission coefficient 

In Equation (13) of their paper, M09 wrote the emission coeffi- 
cient for dust in an optically thick region as: 



T]y = Ky By(T) 



r*oo 

I Ky 

riv — 
Jo Xv_ 

r*co 

rjyd 

Jo 



dv 



CTy 

+ — r] Y 

Xv 



(9) 



Their original equation included a dB v (T)/dT term ins t ead of 
B V (T), because they make use of the lBiorkman & Wood! (l200lh 
temperature correction method for determining the dust temper- 
ature. I use By(T) here for consistency, but throughout the deriva- 
tion, B V (T) can be replaced by dB Y {T)/dT with no impact on the 
final result. 

Equation (0 is similar to Equation (0), but includes an extra 
term which is the ratio of two integrals. It is not clear why this 
term was included by M09, because other than possibly chang- 
ing the dust temperature - which would affect B V (T) - the addi- 
tion of scattering should not affect the thermal emissivity if k v is 
held constant. However, even with this extra term, Equation © 
can in fact be simplified. One can set 



7]y dV 

o Xv 



Jr»oo 
o 



(10) 



dv 



which is a frequency-independent constant for any given dust 
type. Equation (|9]l then becomes: 



l]y = Ky By(T) £+ T]y 

Xv 



M09 suggest solving this through an iterative method, but in fact, 
this can be solved exactly, by re-arranging for rj v as for Equation 
© in 32~fl giving: 



rjy = Xv B Y {T) e 

Substituting this back into Equation (fTOt gives 

Jr»cc 
X v By(T)e — dv 
o Xv 



e - 



\ Xv 

Jo 



(12) 



(13) 



B V (T) e dv 



The e in the integrals cancel out (since they are frequency inde- 
pendent), so that 



e = 



Ky By(T) dV 

_ i^p_ 

> _ 1 

Xp 



f 

Jo 



(14) 



X v B y (T) dv 



where ~kp and xp are the Planck mean mass absorption and ex- 
tinction coefficients respectively. Therefore, Equation (O can 
solved exactly rather than through an interative procedure, with 
the solution given by 



T]V = XV By(T) ^ 

Xp 



(15) 



This is very similar to Equation (01, but includes an extra con- 
stant multiplicative term. However, the derivations in the follow- 
ing sections are valid regardless of whether the emissivity is cal- 
culated using Equation (0 or (fT3T l. as in all cases multiplicative 
constants to the emissivity cancel out. 

2.3. Diffusion coefficient 

The diffusion coefficient is given by M09 as 



D 



(d 2 ) 
6(d) 3 



Jr»oo 
o P 2 X 2 v 



• dv 



Jr>co 
PXv 



(16) 



dv 



where (d) is the mean free path of the photons, (d 2 ) is the mean 
of the path lengths squared, and p is the density of the dust. 
Substituting either Equation © or (Q3) into the above gives 



D 



f 

1 Jo_ 

3 r 

Jo 



' XvByjT) 
n P 2 X 2 y 

' XvB v (T) 



dv 



dv 



(17) 



o PXv 



The Xv terms can be simplified, and the p term can factored out 
as it is not dependent on frequency: 



d --p 



1 Jo 



By(T) 

Xv 



dv 



o 



(18) 



B v (T)dv 



The second term can easily be recognized as 1/xr, the inverse 
of the Rosseland mean mass extinction coefficient. Thus, 



(11) D 



1 



3PXR 



(19) 



This means that the expression for the diffusion coefficient is in 
fact the same whether or not scattering is included. 

2.4. Mean opacity to absorption 

The mean opacity to absorption in the diffusion region is 



Jo 



Ky dV 



o Xv 



f 

Jo 



^dv 



(20) 



Xv 
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This is a mean frequency-dependent mass absorption coefficient 
weighted by the probability of emission at a given frequency 
r] v , and by the mean free path at each frequency, 1 jpx v (where 
the p term cancels out). Substituting Equation (|7]) or ( [13] ) into 
Equation ( l20b gives 



f 



XvB v (T) 

Xv 



Jo 



XvBJT) 

Xv 



(21) 



dv 



As before, this can be simplified by canceling the^- v terms: 
B v (T)K v dv 



f 

Jo 



f 

Jo 



Kp, 



(22) 



B v (T)dv 



which is the standard Planck mean mass absorption coefficient. 



3. Implementation 

In this section, I summarize the M09 MRW algorithm, in the 
light of the new equations derived in ^2] A good criterion for 
deciding whether to start the MRW procedure was suggested by 
M09, and consists in determining whether the distance to the 
closest cell wall is greater than a few times the Rosseland mean 
free path: 



y 

PXR 



(23) 



The y parameter controls the balance of speed versus accuracy. 
If y is set too low, then the sphere used for the MRW may be 
optically thin at long wavelengths, and therefore the photons 
might have escaped directly rather than diffuse if the MRW had 
not been used. This leads to over-estimated temperatures, as de- 
scribed in more detail in M09. If y is set too high, then the per- 
formance gains from using the MRW are decreased, since the 
photons still need to undergo a significant number of interac- 
tions. The starting criterion can be checked after each scattering 
or absorption/re-emission, and does not increase the computing 
time noticeably. 

If the starting criterion is met, a sphere of radius Rq and cen- 
tered on the current photon position can be set up, with Rq being 
at most the distance to the closest grid cell wall. The diffusion 
approximation can then be solved exactly inside the sphere (see 
M09 for details). The diffusion solution leads to the following 
algorithm: first, one samples a random number £ 6 [0, 1] and 
uses it to solve for y in the following equation: 



where D is the diffusion coefficient given by D — 1 /3pxn- The 
energy absorbed by the dust g rains (whi ch is needed to calculate 
the te mperature in both the lLucv|[l999l and iBiorkman & Woodl 
120011 methods) during this MRW is then 



E = Eyd pKp 



(26) 



where E y is the energy of the photon, p is the density, and Kp is 
the Planck mean mass absorption coefficient. 

Finally, the emergent spectrum from an optically thick region 
is given by I v — S v = B V {T), so the frequency of the photons 
exiting the diffusion sphere should be samp led from the Planck 
functio n at the local dust temperature. If the Biorkman & Wood 
(2001) temperature correction method is used, the frequency of 
the photons should be sampled from dB v (T)/dT rather than sim- 
ply B V {T). 



4. Summary 

In this note, I have shown that the MRW equations presented 
by M09 for the emission coefficient of dust, the diffusion co- 
efficient, and the mean opacity to absorption can be simplified 
considerably. The emission coefficient defined by M09 does not 
need to be solved iteratively, but instead is directly given by 
Equation ( fl3l >. The expression for the diffusion coefficient in- 
cluding scattering is in fact identical to that without scatter- 
ing, and is given by Equation ( fT9l . Finally, the mean opacity 
to absorption is simply the Planck mean mass absorption co- 
efficient. All of these values are directly related to the Planck 
and Rosseland mean mass extinction and absorption coefficients, 
which are usually already pre-computed in Monte-Carlo radia- 
tion transfer codes. Thus, the number of calculations involved 
with the MRW can be greatly reduced, and the MRW technique 
can be implemented into existing codes with very little effort. 
An overview of the algorithm, including caveats, is given in |3] 
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<r = 2j](-iry 



(24) 



As y tends to 1, the sum needs to be computed to higher and 
higher values of n in order to preserve a constant numerical ac- 
curacy. The most efficient way to carry out this sampling is to 
pre-compute the sum very accurately for a range of y values, 
and to then interpolate for y given f in the Monte-Carlo code. 
Once the value of y is determined, one can compute the distance 
traveled to exit the diffusion sphere, using 



ct 



(25) 
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